%This program simulates the equilibrium paths leading back to the secular
%stagnation steady state and to the neoclassical steady state.
%Produces Figure 2
tic
clear;
format long;

T=32; %Number of years of horizon

%% Calibration

%Exogenous parameters
rho=0.04; %Discount rate
n=-0.007; %Population growth rate
g=0.013; %Rate of technical progress
frisch=0.85; %Frisch elasticity of labor supply
piW=g; %Downward wage rigidity

%Targets to match
LaborSupply=1; %Labor supply in (neoclassical) steady state
ConsumptionGap=0.1; % -(cSS-cFB)/cFB
MaxPonzi=1.25; %Ratio of Ponzi scheme to GDP in the Ponzi steady state
HalfLife=2; %Years necessary for the inflation anchor to close half the gap with the (constant) rate of inflation
PhillipsSlope=-0.3; %Elasticity of nominal wage inflation with respect to unemployment
MinWealth=2; %(Negative of) minimum wealth relative to consumption under secular stagnation

%Endogenous parameters
kL=LaborSupply^(-1-1/frisch); %Scale parameter of labor supply function
theta=log(2)/HalfLife; %Speed of inflation deanchoring
beta=-PhillipsSlope*frisch; %Slope coefficient of the Phillips curve

%Neoclassical steady state & Secular stagnation steady state
LNeo=LaborSupply;
LSS=(1-ConsumptionGap)*LNeo;
NaturalRate=g-rho*(LNeo/LSS-1)-LNeo/LSS*piW;
lSS=(kL*LSS)^(-frisch); %Labor supply

%Ponzi steady state
deltaPonzi=MaxPonzi*LNeo;

%Calibration of the preference for wealth
RefWealth=-MinWealth*LSS;
sigma=log(LNeo*(rho+piW)/(LSS*(rho-n)))/log(1-deltaPonzi/RefWealth);
kW=((rho-n)/LNeo)*(deltaPonzi-RefWealth)^sigma;

%Taylor rule
pix=0.02; %Inflation target
phi=1.5; %Taylor parameter

%Public debt
S=10; %Maturity of long-term debt
D=1.66; %Long-term indebtedness
H=0.25*LSS; %Helicopter drop

%Initial values
piA0=piW-g;
delta0=H;

%% Solving for the initial conditions of the saddle-path leading to steady state

%Solves for L0 leading to the secular stagnation steady state (through a reverse shooting algorithm starting close to steady state on the saddle-path leading to it)
dd=LSS/(-RefWealth)*sigma*(rho+piW)/(rho+n+2*piW);
pp=0.0001; %distance from steady state
sol=ode23t(@(t,y) dLdDelta(t,y,rho,kW,RefWealth,sigma,piW,n),[pp H],LSS+pp*dd);
%plot(sol.x',sol.y')
L0sim(1)=deval(sol,H);
delta0sim(1)=H;

%Solves for L0 leading to the neoclassical steady state with inflation on target (while solving for the corresponding value of delta0)
Toriginal=T;
L0=1.020939;
dd0=0.069188785340231;
err=1;
while abs(err)>10^(-7)
    delta0=dd0;
    for T=20:Toriginal %To facilitate convergence, the value of L0 is initially solved for T below Toriginal
        SolveL0Ponzi %Solves for L0 such that at time T the inflation anchor piA is on target, equal to pix
    end
    err=H-(1-exp(-interp1(t,y(:,4),S)))*D-dd0; %Gap between guessed value of Delta0 and the value implied by simulated path
    dd0=H-(1-exp(-interp1(t,y(:,4),S)))*D; %Updates the value of Delta0
end
T=Toriginal;
fprintf('Decline in the price of long-term debt (on the path leading to the neoclassical steady state) = %f\n',(1-exp(-interp1(t,y(:,4),S))));
fprintf('Initial magnitude of the Ponzi scheme relative to GDP under stagnation = %f\n',(H-(1-exp(-interp1(t,y(:,4),S)))*D)/LSS);

L0sim(2)=L0;
delta0sim(2)=delta0;

%% Simulating the equilibrium trajectories that are plotted

for k=0:1
    L0=L0sim(k+1);
    delta0=delta0sim(k+1);
    %Solves a differential equation from time 0 to T using ode15s or ode23t
    sol=ode15s(@(t,y) Equi(t,y,rho,g,n,frisch,kL,theta,beta,phi,pix,kW,RefWealth,sigma,NaturalRate,piW),[0 T],[L0;piA0;delta0;0]);
    t=sol.x';
    y=sol.y';

    L(:,k+1)=interp1(t,y(:,1),0:0.01:T)';
    piA(:,k+1)=interp1(t,y(:,2),0:0.01:T)';
    delta(:,k+1)=interp1(t,y(:,3),0:0.01:T)';
    pi(:,k+1)=max(piA(:,k+1)+beta*(kL*L(:,k+1).^(1+1/frisch)-1),piW-g);
    i(:,k+1)=max(NaturalRate+pix+phi*(pi(:,k+1)-pix),0);
end

t=[0:0.01:T]';

%% Figures

FontMain=34;
Font=26;
LW=3.4;
ST=-3;
TT=T-7;

%Consumption
axes('Parent',figure('InvertHardcopy','off','Color',[1 1 1]),'box','on','FontSize',Font,'position',[0.04 0.58 0.43 0.38],'LineWidth',1)
hold on
xlim([ST TT])
ylim([0.88 1.045])
plot([ST;0;t(1:TT*100)],[LSS;LSS;L(1:TT*100,1)],'DisplayName','Secular stagnation','LineWidth',LW,'LineStyle','--','Color',[0 0 0])
plot([ST;0;t(1:TT*100)],[LSS;LSS;L(1:TT*100,2)],'DisplayName','Neoclassical','LineWidth',LW,'LineStyle','-','Color',[0 0 0])
xlabel('Years','FontSize',Font)
title('Output','FontSize',FontMain)
set(legend('show'),'Location','northeast','FontSize',Font,'box','off')

%Inflation
axes('box','on','FontSize',Font,'position',[0.53 0.58 0.43 0.38],'LineWidth',1)
hold on
xlim([ST TT])
ylim([-0.2 2.5])
plot([ST;0;t(1:TT*100)],[100*(piW-g);100*(piW-g);100*piA(1:TT*100,1)],'LineWidth',LW,'LineStyle','--','Color',[0 0 0])
plot([ST;0;t(1:TT*100)],[100*(piW-g);100*(piW-g);100*piA(1:TT*100,2)],'LineWidth',LW,'LineStyle','-','Color',[0 0 0])
%plot(t(1:TT*100),pi(1:TT*100,1),'LineWidth',1.5,'LineStyle','-','Color',[0 0 0])
%plot(t(1:TT*100),pi(1:TT*100,2),'LineWidth',1.5,'LineStyle',':','Color',[0 0 0])
%plot(t(1:TT*100),pi(1:TT*100,3),'LineWidth',1.5,'LineStyle','--','Color',[0 0 0])
xlabel('Years','FontSize',Font)
title('Inflation anchor (%)','FontSize',FontMain)

%Ponzi
axes('box','on','FontSize',Font,'position',[0.04 0.07 0.43 0.38],'LineWidth',1)
hold on
xlim([ST TT])
ylim([-0.025 0.27])
plot([ST;0;t(1:TT*100)],[0;0;delta(1:TT*100,1)],'LineWidth',LW,'LineStyle','--','Color',[0 0 0])
plot([ST;0;t(1:TT*100)],[0;0;delta(1:TT*100,2)],'LineWidth',LW,'LineStyle','-','Color',[0 0 0])
xlabel('Years','FontSize',Font)
title('Ponzi scheme','FontSize',FontMain)

%Nominal interest rate
axes('box','on','FontSize',Font,'position',[0.53 0.07 0.43 0.38],'LineWidth',1)
hold on
xlim([ST TT])
ylim([-0.5 2])
plot([ST;0;t(1:TT*100)],[0;0;100*i(1:TT*100,1)],'LineWidth',LW,'LineStyle','--','Color',[0 0 0])
plot([ST;0;t(1:TT*100)],[0;0;100*i(1:TT*100,2)],'LineWidth',LW,'LineStyle','-','Color',[0 0 0])
xlabel('Years','FontSize',Font)
title('Nominal interest rate (%)','FontSize',FontMain)

toc
